MR Clustering

Notebook for generating mutual rank (MR) clusters

library(boxr)

Welcome to boxr 0.3.4.99999!
Bug reports & Getting Help: https://github.com/brendan-R/boxr/issues

library(magrittr)
library(knitr)
library(limma)
library(snowfall)

Loading required package: snow

library(igraph)

Attaching package: ‘igraph’

The following object is masked from ‘package:magrittr’:

    %>%

The following objects are masked from ‘package:stats’:

    decompose, spectrum

The following object is masked from ‘package:base’:

    union

box_auth()

Reading client id from .Renviron
Reading client secret from .Renviron
Auto-refreshing stale OAuth token.
boxr: Authenticated at box.com as Julin Maloof (jnmaloof@ucdavis.edu)

box_setwd(24503987860)

box.com working directory changed to 'BrapaNetworks'

      id: 24503987860
    tree: All Files/BrapaNetworks
   owner: jnmaloof@ucdavis.edu
contents: 9 files, 0 folders

box_ls()

box.com remote object list (9 objects)

  Summary of first 9:

                            name type           id   size
1 Brapa2011BayesHeight_blups.csv file 164446516413  15 kB
2 Brapa2012BayesHeight_blups.csv file 164446639332  28 kB
3    Brapa_V1.5_annotated.csv.gz file 163054185089 2.7 MB
4   brassica_gene_clusters.RData file 164456845216 110 kB
5               CountsVoom.Rdata file 162362592785 240 MB
6         NormalizedCounts.RData file 162341539458 100 MB
7                RIL_TMM_CPM.csv file 162352229544 630 MB
8           RIL_v1.5_mapping.tsv file 162218089685 110 MB
9                     test.Rdata file 162921069376   64 B
                 owner
1 jnmaloof@ucdavis.edu
2 jnmaloof@ucdavis.edu
3 jnmaloof@ucdavis.edu
4 jnmaloof@ucdavis.edu
5 jnmaloof@ucdavis.edu
6 jnmaloof@ucdavis.edu
7 jnmaloof@ucdavis.edu
8 jnmaloof@ucdavis.edu
9 jnmaloof@ucdavis.edu


Use as.data.frame() to extract full results.

First attempt, use RIL means

Get data

annotation  <- read_csv("~/Box Sync/BrapaNetworks/Brapa_V1.5_annotated.csv.gz")

Error: could not find function "read_csv"

Summarize means

save(voom.fit,fil="~/Box Sync/BrapaNetworks/voom.fit.Rdata")

Error in save(voom.fit, fil = "~/Box Sync/BrapaNetworks/voom.fit.Rdata") : 
  object ‘~/Box Sync/BrapaNetworks/voom.fit.Rdata’ not found

RIL.coefs <- grep("RIL",colnames(coef(voom.fit)))
RIL.pvals <- topTable(voom.fit,coef = RIL.coefs,number = Inf)
sum(RIL.pvals$adj.P.Val<0.05)

[1] 26721

sum(RIL.pvals$adj.P.Val<10e-10) 

[1] 16596

So pretty much everything is highly significant

Start by taking the top 10000

genes.of.interest <- row.names(RIL.pvals)[1:10000]
expr.data <- voom.fit$coefficients[genes.of.interest,]
expr.data <- expr.data[,-grep("trt",colnames(expr.data))]
dim(expr.data)

[1] 10000   124

head(expr.data[,1:6])

          (Intercept) pdata$RILRIL_103 pdata$RILRIL_104 pdata$RILRIL_113
Bra001964    1.082235       -4.0321896     -4.091759090        3.6994071
Bra002728    5.767475       -0.5505380     -0.380674984       -3.2039164
Bra002785    7.202490       -0.4198373     -0.037837754       -1.9963410
Bra004960    4.733936       -3.8185006     -0.004181862       -0.2564264
Bra006384    4.969959       -6.8228958     -7.276735466       -0.3784078
Bra006411   -2.499584        6.3062481     -0.470837158        0.1143828
          pdata$RILRIL_115 pdata$RILRIL_12
Bra001964       0.01722016      0.64266465
Bra002728      -0.60526951     -0.81853556
Bra002785      -0.04677711     -0.05993573
Bra004960       0.14101788     -0.22825526
Bra006384      -8.39926386     -7.41766487
Bra006411      -0.72175621     -1.10891611

expr.data[,-1] <- expr.data[,-1] + expr.data[,1] # add the intercept to the coefficients
colnames(expr.data) <- sub("pdata$RIL","",colnames(expr.data),fixed=TRUE)
head(expr.data[,1:6])

          (Intercept)    RIL_103   RIL_104   RIL_113   RIL_115    RIL_12
Bra001964    1.082235 -2.9499545 -3.009524  4.781642  1.099455  1.724900
Bra002728    5.767475  5.2169367  5.386800  2.563558  5.162205  4.948939
Bra002785    7.202490  6.7826522  7.164652  5.206149  7.155712  7.142554
Bra004960    4.733936  0.9154359  4.729755  4.477510  4.874954  4.505681
Bra006384    4.969959 -1.8529365 -2.306776  4.591552 -3.429305 -2.447706
Bra006411   -2.499584  3.8066642 -2.970421 -2.385201 -3.221340 -3.608500

Get Rob’s Data

blups2011 <- box_read_csv(164446516413)

Remote file '/var/folders/xr/9cbydt955pj42zfq6mc_y5g40000gn/T//Rtmp7o2lKa/Brapa2011BayesHeight_blups.csv' read into memory as an object of class data.frame

blups2012 <- box_read_csv(164446639332)

Remote file '/var/folders/xr/9cbydt955pj42zfq6mc_y5g40000gn/T//Rtmp7o2lKa/Brapa2012BayesHeight_blups.csv' read into memory as an object of class data.frame

blups2011$Line %<>% paste("RIL",.,sep="_")
blups2012$Line %<>% paste("RIL",.,sep="_")
blups2011 %<>% subset(select=!grepl("V1|individual|blk|trt|treat",colnames(.)))
blups2012 %<>% subset(select=!grepl("V1|individual|blk|trt|treat",colnames(.)))
head(blups2011)

head(blups2012)

Merge Data

expr.blup.2011 <- merge(blups2011,t(expr.data),by.x="Line",by.y=0)
expr.blup.2012 <- merge(blups2012,t(expr.data),by.x="Line",by.y=0)
head(expr.blup.2011)

head(expr.blup.2012)

Correlation

expr.blup.2011.cor <- cor(expr.blup.2011[,-1])
expr.blup.2012.cor <- cor(expr.blup.2012[,-1])

Mutual Rank

sfStop()

Stopping cluster

save(blups2011,blups2012,voom.fit,expr.blup.2011, expr.blup.2012,expr.blup.2011.mr, expr.blup.2012.mr,file = "~/Box Sync/BrapaNetworks/MR.Rdata")

Get graphs at different MR

First: define some functions

Get blup networks at different MR thresholds

Use different MR thresholds of 10 to 210 (absolute values) to build the network (genes with a correlation higher than the threshold were considered connected). I then subset the network to take all genes “directly” connected to LFY.

plots

null device 
          1 

null device 
          1